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We analyse the dynamics of a hard-sphere lattice gas on generalised Bethe lattices using a projec- 
tive approximation scheme (PAS). The latter consists in mapping the system's dynamics to a finite 
set of global observables, closure of the resulting equations is obtained by approximating the true 
non-equilibrium state by a pseudo-equilibrium based only on the value of the observables under con- 
sideration. We study the liquid-crystal as well as the liquid-spin-glass transitions, special attention 
is given to the prediction of equilibration times and their divergence close to the phase transitions. 
Analytical results are corroborated by Monte-Carlo simulations. 

PACS numbers: 



I. INTRODUCTION 



The statistical mechanics of finite-connectivity systems has seen a large increase in interest within the last decade. 
The major reason for this interest was the emergence of a close interdisciplinary collaboration between statistical 
mechanics and theoretical computer science. In many hard combinatorial problems, phase transitions in resolvability 
and algorithmic performance were observed - once these problems are suitably randomised Q,Q>IS- This observation 
obviously provides a large temptation to use statistical mechanics tools to get insight into the problem structure going 
beyond the one which can be obtained with traditional tools from discrete mathematics and computer science. The 
natural mapping of these models to disordered finite-connectivity models was therefore used to investigate problems 
like satisfiability 0,03, El; vertex cover 0, graph colouring Q etc. 

Coming back to physical models, finite-connectivity systems can be understood as an intermediate step between 
fully-connected mean-field models |9j and realistic, finite-dimensional systems. Even if their structure is not geometri- 
cal in the sense that they are embeddable into any finite-dimensional space, the finite connectivity allows for a reason- 
able implementation of notations like neighbourhood, distance etc. This property enables construction and analysis of 
microscopically motivated models, in particular for spin-glasses 0,0] and structural glasses [l3l IT3 . [TEl fla . fl7 | . 

The understanding of these models is, however, mainly based on the analysis of their thermodynamic equilibrium 
behaviour. Due to a technical break-through based on an application of the cavity method to finite-connectivity 
models 0], the structure of stable and metastable states, the appearance of phase transitions etc. are relatively 
well-understood by now. 

The knowledge is much less evolved in what concerns the dynamics of finite-connectivity models [Til flti, I20I I2TL I22I 
l23l | . Here we will use an approximation scheme first introduced to analyse the dynamical behaviour of stochastic local- 
search algorithms for combinatorial optimisation |24ll25j |. We will, however, use a different formulation which is based 
on the observation from Ref. |22| that the method can be put into the canonical form of a projective approximations 
scheme being equivalent to the dynamical replica theory [2(J, [27], [28j which was originally developed for fully connected 
systems. Note also 29], where the same method was recently applied to Ising models on random graphs. 

The model studied here is a hard-sphere lattice-gas model on a generalised Bethe lattice, which was first proposed in 
|14| as a microscopically motivated, but solvable model for the structural glass transition. The equilibrium behaviour 
of this model, including crystallisation as well as glass formation, was studied in great detail in 0]. The present 
work, which is concentrated on the relaxational dynamics of the model, can be seen as a natural continuation of the 
work presented therein. Even if we have tried to make the present paper as self-contained as possible, sometimes we 
have to refer to technical details in the previous publication - whose full inclusion would go beyond the scope of the 
present article. 

The paper is organised as follows. Sec. ITT1 gives an introduction to the system including a brief review of its static 
properties and the definition of the dynamics. In the following section, the general outline of the applied approximation 
scheme is presented. Sec. IV is dedicated to crystallisation dynamics, containing two approximations and comparison 
with Monte-Carlo simulations. In Sec. V an approximation of the dynamics near the spin-glass transition is presented. 
Finally, conclusion and outlook are given in the last section. 



FIG. 1: Part of a generalised Bethe-lattice with k = 2, p = 3. The 21 cliques connecting the central vertex to its nearest and 
second neighbours are shown. Figure from Ref. |l7j|. 



II. THE MODEL 

In a first step, we are going to define the model itself, review its static properties which were analysed in detail in 
|17| . In the second part of this section we will define microscopic physical dynamics for the model, whose analysis will 
be the subject of the whole paper. 

A. Construction and static properties 

The system is a lattice gas of hard spheres living on a generalised Bethe lattice. Let us first specify the character 
of the particles. The property of hard spheres is realised by introducing an excluded volume associated with each 
particle. This is done by requiring that two neighbouring sites (i.e. sites that are connected with an edge in the graph) 
cannot be occupied simultaneously. Two particles have thus a minimal distance of two on the lattice. 

The generalised Bethe lattice is constructed in the following way: The basic building blocks are small completely 
connected sub-graphs (i.e. cliques), each of which contains p + 1 vertices. In each vertex k + 1 of these cliques are 
merged. This procedure leads to a graph which is locally isomorphic to that shown in Fig. ^ for ° ne specific choice 
of k and p. The ideal graph can be thought of as an infinite continuation of the shown structure such that the whole 
graph remains loop-less apart from the inner-clique loops. Note that the ordinary Bethe lattice can be obtained as the 
special case p = 1. It is furtheron obvious that a close packing on such a lattice consists of a configuration where each 
clique carries exactly one particle - which, because of the hard-sphere constraint, realises also the maximum particle 
density that can be found in a single clique. 

A treatment of a lattice gas on this graph would, however, lead to subtle questions about the boundary conditions. 
In the thermodynamic limit, we would face a strange situation where the boundary is infinitely far from a fixed centre, 
but the number of boundary sites constitute a finite fraction of all sites. As we shall recur to finite lattices when 
confronting our results with Monte-Carlo data, it seems more appropriate to base the calculations on the properties 
of finite random realisations of the generalised Bethe lattice from the beginning, see also the discussion in [l7| . 

A finite graph allowing for a close packing as described above can be constructed, e.g., as a bipartite (or (p + 1)- 
partite) graph which possesses one (or p+l) close packings. Any given close packing defines two sub- lattices referred 
to as 0-lattice (empty sites) and 1-lattice (occupied sites) in the following. The equilibrium properties of the system 
can be obtained from a calculation of the grand canonical partition function 

E{fx)= J2 e^™^ Yl (l-n inj ) (I) 

ni,..,n»e{0,l} {i,j}eE 

where /z is the chemical potential, i, j = 1, . . . , N label the sites and E is the set of edges of the graph. The calculation 
was carried out in detail in 17]. One finds that the system has a homogeneous low-/i-phase (i.e. both sublattices have 
same particle densities) which is referred to as the liquid phase. In the high- /i- regime the system can be found in 
a phase where the 1-lattice has a higher particle density than the 0-lattice. This phase is called crystalline. In the 
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case p > 1 there is also an intermediate regime where both liquid and crystalline phases are locally stable, and the 
phase transition between liquid and crystal is of first order. This explains the interest of generalising the concept of 
an ordinary Bethe-lattice (p — 1) which lacks the intermediate phase. The extension can be motivated physically by 
the abundance of short loops found in finite dimensional lattices. They are taken into account at least approximately 
by the introduction of local cliques. 

Another possibility to realise the finite graph is to construct a generalised regular random graph, i.e. a graph that 
is locally isomorphic to a generalised Bethe lattice but contains long loops with a typical length that is logarithmic 
in N. This setting leads to a qualitatively different high-/i-phase: The graph is, with high probability, no longer 
compatible with a cry stalline close packing. The statics of the system can be appropriately calculated via the cavity 
approach, cf. [Tol Il6| . The scenario for arbitrary k and p is described in fill . In the present study we only focus 
on the case where p — 1 (ordinary Bethe lattice). The system then displays a liquid low-/i-phase and a continuous 
transition to an amorphous high-// spin-glass phase. A one-step-replica-symmetry-broken (1RSB) ansatz takes into 
account the frustration due to odd loops, and gives a lower grand-potential than the liquid solution. 



B. Definition of the dynamics 

The dynamics of the model is made up by two basic processes: The first one being diffusion of particles on the 
lattice, the second one particle exchange with an infinite particle bath. The realisation of these processes is as follows: 

• In every time step, a lattice site is selected randomly. 

• If the lattice site is occupied, with probability qo the particle is annihilated (or, equivalently, transferred to the 
particle bath). With probability q s , the particle tries to jump to a randomly selected neighbouring lattice site. 
This jump is only realisable, if the new lattice site is empty itself, and has no further occupied neighbouring 
sites. Every other case would either result in two particles on one site, or in two neighbouring particles, which 
is excluded according to Eq. Q). 



If the lattice site chosen in the first step is empty, we try to create a new particle with probability qi . This trial 
is successful only if also all neighbouring sites of the selected one are empty. 



The Monte-Carlo (MC) dynamics defined in this way approaches the thermodynamic equilibrium described by Eq. 
if the rates for particle annihilation and particle creation are related by detailed balance: 

ef = 21 (2) 

which we assume to be satisfied throughout this paper. Besides the fact that the described dynamics can be easily 
implemented on arbitrary finite lattices or graphs, it can also be understood as a valid microscopic dynamics of the 
physical model under consideration. 

The time unit is one MC sweep (At = 1) which corresponds to a number N of the above trials. A single trial thus 
corresponds to At = 1/N, and the time variable becomes continuous in the thermodynamic limit. 



III. PROJECTIVE APPROXIMATION SCHEMES 

The dynamics of diluted systems have been described successfully in the past by use of rate equations [22, |3 EH EH- 
This method, even if formulated differently in some of these references, has an elegant formalisation as a projective 
approximation scheme (PAS) described in the following steps: 

• A set of intensive macroscopic observables O is chosen. Typically O are densities and we require "O 3 O equ u v , 
i.e. the set must contain the observables necessary for the discrimination of the different phases at equilibrium 
(e.g. the two sublattice densities to describe crystallisation in our model). 

• In a general non-equilibrium situation, the dynamical equations for the observable do not close in the set O. 
One can, however, approximately close the equations by assuming equilibrium in a generalised ensemble: All 
microscopic configurations leading to the same values of all observables in O are assumed to be equiprobable. 

There is, of course, no reason why this closure assumption should be true in general. It can, however, be seen as the 
most natural approximation to the actual distribution of the system which is based solely on the knowledge of the 
observable values, and on no further information. In this sense, the PAS has to be considered a the best possible 
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approximation of the model's dynamics based on O. This does not mean that, for specific initial conditions, parameter 
values etc., there cannot be a better approximation to the actual time dependence of the observables. For the general 
case, it constitutes, however, an optimal use of the available information on the system. Further on, since we required 
u O 2 Oequu", we are sure that the true thermal equilibrium is a fixed point of the PAS. 
Technically, the equipartition hypothesis is realised in the following way: 

• For each observable in 0, a conjugate chemical potential is introduced (e.g. the usual chemical potential p 
would couple to the particle density p). The potentials are collected in the vector p. 

• The generalised grand partition function 

3 G (/Z) = 5>"* d J] (l-mnj) 
{«;} {i,j}eE 

automatically implements the equipartition condition. 

• For given observable values of O, the values of the chemical potentials p are obtained self-consistently from 

s g (m). 

• Starting from the values of the observables O, we derive equations for the time changes O. The latter depend in 
general on a larger set of observables, whose values are evaluated in the generalised equilibrium ensemble using 
the previously determined chemical potentials. In the next chapter it will become clear, how this step is realised 
technically 

This scheme provides a well-defined prescription for the approximation of the dynamics of the system based on the 
observable set O. Unfortunately this prescription does not include any intrinsic way of quantifying deviations from 
equipartition, such that all results based on this scheme have to be cross-checked against numerical simulations. On the 
other hand, it is obvious that going to more detailed observables means to enhance the quality of the approximation, 
at the cost of introducing a more complex generalised ensemble. This will be shown in detail for the crystallisation 
dynamics below in Sec. II VI 

We mention that the assumption of equipartition of probability was introduced in the context of dynamical replica 
theory [2g for fully connected spin-glass models. In this context, the technical treatment of B.g(0) was obtained using 
the replica trick. In the context of finite-connectivity systems, however, the application of the cavity method will be 
much more efficient. 



IV. CRYSTALLISATION DYNAMICS 

As a first application of PAS to our model, we are going to study the crystallisation dynamics of the model. We 
will show two levels of PAS, the first one being the minimal description leading to a correct characterisation of the 
equilibrium, the second one will include a larger set of observables. Both schemes will be compared to numerical data 
obtained from MC simulations on large, but finite lattices. The experience gained in this context will be used in a 
later section in order to approximate also the dynamics close to the spin-glass transition. 

A. The p-approximation 

From the generalised ensemble to approximate dynamical equations 

Let us start with a minimal set of observables. As we need an indication whether the system is liquid or crystalline, 
we obviously need the two sublattice particle densities p' ' and p^ that refer to the 0-lattice and 1-lattice respectively. 
The normalisation is chosen such that p^ — ► and p^ — > 1 in the close packing limit p, — ► oo. Before being able to 
write any closed dynamical equation for the two sublattice densities, we have to consider the generalised equilibrium 
ensemble defined in the last section. We introduce two chemical potentials p^ and p^ that couple to the particles on 
the respective sublattice. The calculation of Sc can be performed analogously to that of the grand partition function 
using the concept of rooted trees. Here we give only a shortened description of how to do so, for technical details 
see [13. A rooted tree is obtained from a generalised Bethe lattice by choosing a site as the root and removing one 
of the k + 1 incident branches. Rooted trees can be iterated according to Fig. [3 where dashed lines depict the edges 



FIG. 2: Possible iteration steps leading to rooted trees with root in the 0-lattice (left), 1-lattice (right). Vertices belonging to 
the 0-lattice (1-lattice) are depicted by o (□). Figure from Ref. [r^|. 



which are added in one iteration step. We introduce conditional partition functions (resp. sl°^) for a subtree 

rooted in a 0/1-site which is fixed to be empty (resp. occupied). These quantities iterate according to 

s£°>' = (ffa^^^awjf^ + ip-i)^^)* (3i) 

S (0)' = j» (zW(E?)y-i) k (3ii) 
S^)' = ((E^Y+pE^iE^y-'f (3iii) 
- e^((E^y) k . (3iv) 



The stroke labels the partition functions of the resulting large rooted tree. All 0-rooted trees that enter in the 
iteration are taken to be identical, the same is assumed for all 1-rooted trees. When we pass to intensive quantities 

by introducing local fields cc' / 1 ) = ^jom > the strokes in Eqs. © drop out and the iteration equations reduce to two 

equations for the local fields x^ / 1 - 1 which still contain the unknown chemical potentials p^ / 1 ^. These can be eliminated 
through the connection of Eqs. 10 with the observables p 1 ^/ 1 ^ which is established via the physical partition functions 
si / 1 ) and 3* . The latter refer to a site within the generalised Bethe lattice (not the root of a subtree), i.e. they 
are obtained like the stroked partition functions in Eqs. Q after the substitution k — ► k + 1 on the r.h.s. of the 

(o/i) S(o/i) 

equations. The relations jtrwr; = ^*o/i) allow to derive two more equations, and we end up with four equations 
containing the four unknown a;' / 1 ) and p(°^\ Their solution is straight-forward yielding 

,(0/1 
l-p~ pp (0) 

and 

(0/1) p {o/i)(l- pW))* 



x (o/i) P_ (4) 



(5) 



We note that Eq. JSJ is the self-consistency result for the chemical potentials and Eq.0|is the result for the conditional 
partition functions that are essentially contained in the local fields. We also mention that by setting p^ = fiyl = p 
the results of the calculation for the statics are recovered, i.e. we have found an "equilibrium" treatment for non- 
equilibrium configurations of p^ ' 1 ' that naturally contains the physical equilibrium ensemble as a special case. 

In order to derive the approximative differential equations for the time-evolution of p^ / 1 ), we start with the exact 
equations for the average change AN^^ 1 ^ in the number of particles A^i ^ 1 ' on the respective sublattice. The average 
is taken over the time At = l/N where N = + is the overall number of sites and N^- ^ 1 ^ are the numbers of 
sites belonging to the sublattices. At is the time that corresponds to one action of the MC-algorithm that implements 
the dynamics (see Sec. Ill The changes AAi ^ are given by rate-equations: 

aJ V<°/« = ,<»A) . + ^(1 -/>/») ,M" m 

- V- pi0/1) q * ^ 0/1) + V p{1/0) qs v ' im ■ (6) 
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The first term on the r.h.s. corresponds to a loss due to annihilation of a particle. The three factors mean that the 

(0/1 
W 



algorithm must choose a site on the lattice in question ( Ni J - ) which is occupied (/^ 0//1 ^) and annihilate it (qo). This 



contribution is exactly known. The second term contains the gain due to particle creation. Here •p^ ^ denotes the 
probability that the particle creation is possible without violating the hard-sphere constraint, i.e. the probability that 
all neighbours of the selected site on the 0/1-lattice are also empty. Its value in one microscopic (non-equilibrium) 
configuration cannot be expressed as a function of the //"/^-values of this configuration, as it contains information 
about nearest-neighbour correlations that are not included in the single-site observables p^^K The third term gives 
the loss of particles caused by those that jump to the other sublattice whereas the last term contains the gain of 
particles caused by those that jump from the other sublattice to the sublattice under consideration. Hereby, 7 3 i° / ' 1 ' > 
denotes the probability that a tried move of a particle on the 0/1-lattice to a neighbouring site on the 1/0-lattice is 
in fact feasible. This is again a quantity which involves neighbour correlations and cannot be expressed exactly for 
one configuration given only by p^ / 1 ) . 

In order to close Eq. JBJ we now proceed to the last step mentioned in the outline of PAS. We have to calculate the 
probabilities 7^°^ within the pseudo-equilibrium ensemble characterised by the already determined formal chemical 
potentials p^°^ l \ i.e. the T^ / 1 ) are calculated via a flat average over all configurations having the desired values of 
pi / 1 ). For instance, we can calculate 

all neighbouring sites are empty 



h-\- 1 

,(0) 



( H « (s (o ))P -i + S (D (s (o ))P -i + {p _ 1) S (i) s (0) (S (o )r - 2 ) 

V v 

no requirements to the neighbouring sites 

l 



fe+i 



(l + + {p - 1) xW) k+1 

® fl-pw-ppl^ k ' 



1 - p(°) 



(7) 



Similarly we find 



For the calculation of we must consider that, starting from a 0-site, only k + 1 of the (k + l)p neighbouring 
sites belong to the 1-lattice which is taken into account with a factor i. Furthermore we know that the ending site 
of the jump is empty as are all other sites in the clique that is shared by starting site and ending site. The only 
condition for is that the kp sites of the k other cliques containing the ending site are empty. These sites are all 
in the 0-lattice as the ending site is in the 1-lattice. This leads to 

k 



V (0) = I 



(N 0) r) _i (i-p^-p^\ k 



p ((E^r+p-^iE^y-y l > 



(9) 



and similarly 



1 - 

These approximations close Eq. ©. We add that 



A- 



a ly(QA) _ ^ (0/1) A(^ (0/1) /iV (0/1) ) _ jV (0/1) Apt"/ 1 ) 

* ~ N At ~ N At [ ' 



and 

N (o) p N (i) ! 



(12) 

N p+1 N p+ 1 v ' 
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for the generalised Bethe lattice. 

Using these ingredients, in the thermodynamic limit N — ► oo Eq. © becomes: 

^=-p<°>» + (1-^)91 ^"^"^ 



1 - p(°) 



,<« ifi^z^, (13i) 



and 



- p n i-p(o) J i-pd) J • (l3u) 

Eqs. constitute what we call the p-approximation, i.e. the approximation for the crystallisation dynamics which 
is derived for a minimal set of observables - p^ and pW - i n the context of PAS. The approximative time-flow of 
pi / 1 ) can be obtained for given initial conditions by integrating Eqs. numerically. Results are shown in Sec. II V 01 
where the p-approximation and <7j-approximation are compared with MC simulations. 



Asymptotics of the p- approximation 

Exactly at the stationary point of these equations, corresponding to an equilibrium solution of the model, the 
assumption of the PAS becomes true. The stationary points of Eqs. are therefore identical to the static liquid 
and crystalline solutions discussed in Ref. [l7j . 

A linear expansion of Eqs. therefore allows to analyse the stability of the equilibrium solution as well as to 
determine the equilibration time. To achieve this, we have to calculate the Jacobian of the above system of equations, 



This matrix has to be evaluated at the stationary point. Let us assume that it has the two eigenvalues Ai and A2 
with corresponding eigenvectors v\ and V2 ■ The time evolution of a small perturbation around the stationary point is 
then given by 



p(°>(t) 




ai e Xlt vi + a 2 e X2t v 2 (15) 



with constants ai and ct2 which follow from the initial condition, i.e. from the specific form of the perturbation. 

The stationary point is obviously stable if both eigenvalues Ai and A2 are negative. In this case, the equilibration 
time is given by r eq = — 1/ max{Ai, A2}. If at least one of the eigenvectors is positive, the stationary point is unstable 
and thus physically irrelevant: Any small perturbation having a component in the direction of the corresponding 
eigenvector carries the system exponentially fast away from the stationary point. 

Let us first consider the liquid solution, which is characterised by equal sublattice densities 

p (0) = p (i) = p 

In this 

case, the eigenvalues can be given in a compact form: 

Ai = - fta + rffr- 1 )) . (16i ) 



ef(l-(p + l)p)(l-p) 

: + l)p) (qiP + p(p 
e*»p(l-(p+l)p)(l-p) 



^ = (1 - (fc + l)p) (gip + p(p + l)(g. - q lP )) 



and the eigenvectors are 



vi=[{), (17) 
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Both eigenvectors are negative for low densities, but A2 changes sign at p cr — This point coincides with the local 
instability of the liquid state in the statical calculation of |l7j ■ 

It is also interesting to look at the equilibration time r eq = — l/max{Ai, A2}. For small chemical potential p, the 
latter is given by — 1/Ai. Looking at the eigenvector vi, we see that the corresponding slowest equilibration process 
in the system describes the decay of global density fluctuations. Close to the instability point of the liquid solution, 
we find, however, r eq — — I/A2, and the eigenvector indicates that the slowest decaying perturbation is a difference in 
the two sublattice densities, i.e. a fluctuation toward a crystalline state. At a certain p = p eq the two eigenvectors are 
equal. At this point, a crossover between the two different slowest processes takes place. It is also the point, where 
the equilibration time takes its minimum inside the liquid phase, i.e. where the system equilibrates faster than for any 
other p < fi cr , see also the discussion in Sec. II V 01 where these analytical results are compared with MC simulations. 

Let us now discuss the stability of the two crystalline solutions calculated in [l^ ■ The determinant of the Jacobian 
matrix A vanishes at the point where the crystalline solution appears. 

• For p = 1, i.e. on the ordinary Bethe lattice, the two crystalline solutions appear continuously out of the liquid 
solution, and they become both stable, marking thus a second-order transition. The two solutions are in fact 
identical, they are related by an exchange of the two sub-lattices which are isomorphic for p = 1. 

• For p > 1, the two crystalline solutions appear discontinuously. The normal crystalline solution, i.e. the one with 
higher total density, becomes immediately locally stable, both eigenvalues of the Jacobian matrix are negative. 
The second solution is first locally unstable and thus unphysical. At p = fj, cr , the two sublattice densities cross, 
however, at the value p cr — an d thus coincide with the liquid solution. For higher p we have po > pi, in 
|17| we have introduced the name "inverse crystallisation" for this phenomenon. This solution takes over the 
local stability of the liquid solution as can be inferred from the eigenvalues of A. 

It is also possible to analytically characterise the critical behaviour close to p cr , which is, on both sides of the 
transition, characterised by a critical exponent — 1: 

T eq <x\p- Pcrl^ 1 ■ (18) 

We skip the calculation of the prefactors here, we only mention that in the case p = 1 both sides of the transition 
are related by an universal factor 1/2, as already found in the case of a simple ferromagnet on a Bethe lattice [22||. 
Further results are given below in the context of the comparison between analytical and numerical data. 



B. The ctj -approximation 

The p-approximation was derived according to PAS using the pure single-site observables p' / 1 ). There we have 
seen that it was necessary to approximate nearest-neighbour correlations within the generalised ensemble. We now go 
over to more detailed observables including also such nearest-neighbour correlations. This can be achieved by defining 



o~j := density of sites on the 0/1-lattice which are empty and have ^g) 
exactly j occupied neighbouring sites, j = 0, . . . , k + 1. 

Note that the cj° imply the previous observables via p^ / 1 ) = 1 — Ylj=o a j ■ Therefore they allow for an exact 
characterisation of the liquid and crystalline equilibrium states, and according to the rule of our PAS they represent 
a valid set of observables. 

Again, we apply PAS introducing formal chemical potentials fj,^ 1 ^, j = 0, . . . , k + 1, being conjugated to the crj ^ 1 '. 
This allows to calculate conditional partition functions aj ' , j = 0, . . . , k which stand for a 0/1-rooted tree with an 
empty root having j occupied neighbouring sites, and S* which stands for a 0/1-rooted tree with occupied root. 
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The iteration of rooted trees (see Fig.^J reads in terms of the conditional partition functions: 




S (°)' = (f^^e-^S^fs^e-^^)" 1 } , (20i) 



1=0 



+ (P 1) Sl 0) ( £ ^.e-ft W f) f £ e^e^'sM 



2=0 7 M=0 




K 2=0 



fe \ p—1 



1=0 



"(1)' 



.7 





k-j 



and 

B W' - J I^A"^-^) I , (20iii) 



;){^(E^.^r} J x{(^)-p. 



The stroked partition functions belong to the rooted tree that results from the iteration. We shall only explain the 
structure of Eq. (|20ii(l exemplarily. It describes the situation where the root of the resulting tree is in the 0-lattice (see 
Fig. |21 left) and empty with j occupied neighbouring sites. The exponential prefactor accounts for the addition of a 
site of this kind. The binomial factor counts the possibilities to distribute the j particles on the k adjacent branches. 
The first curly brackets refer to the situation of a branch containing a particle. This particle can sit on the 1-lattice 

(first term in the brackets, H^). The factor following si 1 -* expresses the possibilities for empty sites with a given 

(o) _ (o) 

number of occupied neighbours on the other p—1 sites of the clique. Note that the factor e t ^ l + 1 e is due to the 
particle present on the 1-site of the clique which increases the number of occupied neighbours for the empty 0-sites 
by one. The second term in the brackets stands for the p—1 possibilities of placing the particle on a 0-site. The 
second curly brackets refer to an empty clique. All possible combinations of empty sites with I = 0, . . . , k occupied 
neighbouring sites are included. 

For the calculation of the partition functions we follow the same lines as for Eqs. ©. We introduce intensive 
quantities with the local fields 

CT (o/i) 

4 0/1) = 3o7i)' *=<),...,*. (2i) 

I — * 

By taking the ratio of Eq. H20ii|) and Eq. (|20i[l as well as Eq. (|20iv() and Eq. I|20iii[) for j = 0, . . . , k we obtain 2(fc + 1) 
equations for these local fields that still contain the 2(k + 2) chemical potentials. These can be eliminated via the 
introduction of physical partition functions sj ^ 1 -', j = 0, . . . , k + 1 and sl ^ which replace the stroked quantities in 
Eqs. H20f> after substitution of k — ► k + 1 in the binomials and exponents. By taking the ratios of Eqs. (|20|l as before 

§(0/1) a (0/l) 

but for j = 0, . . . , k + 1 and using that T| 0/1) = ( /i) we obtain 2(fc + 2) more equations, so that we end up with the 

same number of equations and unknown local fields and chemical potentials. 

The solution of the corresponding equations for the p-approximation was obtained by simple algebra which seems 
not possible in the present case considering the complexity of the equations. We can, however, show that the ansatz 

(o/i) = a( o/i)fc+iii (o/D (22i) 
J k + 1 J 

e-»f™ xfl*> = ^ J -±±a^P (22ii) 
k ~\~ 1 



with j = 0, . . . , k solves the equations under certain conditions. One of coefficients c/ / 1 ) and /J*- / 1 - 1 can be chosen 
arbitrarily for each sublattice. These two extra degrees of freedom arise from the consistency conditions that must be 
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imposed for the validity of Ansatz (|22[1 and which read 

k t , -. k j i -i -, k 



l £+1 (i) (i) _ i±l (o) _ P^l j±i (i) f23l 



z=o 2=0 ;=o 



and which reduce the independent observables by one for each sublattice. Eqs. (|23|l can be easily understood through 
the observation that, e.g., empty sites on the 1-lattice with occupied neighbouring sites contribute to the particle 
density on the 0-lattice (first equation). Eqs. ll'-il) are indeed important conditions as they assure that a choice of 
crj / 1 -* is realisable within a microscopic configuration, and thus they do not imply any restriction of Ansatz (1221) for 
physically relevant situations. 

As in Sec. IIV Al we now derive the differential equations for the time evolution of the observables. Unknown 
probabilities are again approximated by their pseudo-equilibrium values in the ensemble characterised by Eg- A 
complete list of the contributions due to the different actions of the dynamics is given in Appendix ^ We shall only 
discuss the case of the jump from a 0-site to a 1-site which qualitatively includes the calculations occurring in all other 
cases. 

The occurrence of a 0-1-jump requires that the algorithm selects a 0-site (1) which is occupied (2) and the decision 
to try a jump must be made (3) in favour of an ending site which is in the 1-lattice (4). The probability is given by 
(1) • (2) • (3) • (4) = • ■ q s ■ |. ft remains to calculate the probability that in the above situation the jump 

will be successful. To this purpose we introduce the quantity E'} ^ 1 ^ = Ya=q e - ' 1 ' 1 'si / 1 ' which corresponds 

to a rooted tree with an empty root but with the modification that one extra occupied site is connected to the root 
(empty sites connected to the root will not alter the expression). Then the situation where (1) to (4) holds can be 
described by 

3 (o) (s /(o) )P -i 3 / (1)) (24) 

in terms of conditional partition functions. The situation is depicted in Fig. [3] where for grey sites it is not known 
whether they are empty or occupied. For a successful jump we must require all sites labelled with f to be empty. This 
is a restriction of the above situation which reads in terms of conditioned partition functions: 

sffSff-'e^'e-^f . (25) 
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The ratio of Eqs. 125|) and (|24[) gives the probability that the jump can actually be performed: 

Si 0) (E^Y- 1 exPe-Wsg) _ e^e'^E^ 



S (0) (5 ;(0) )P _ 1S /(D Eto^e-"< (1) 3« 

u (1) -u (1) (1) , , 1 



(26) 



i+i 



We shall now discuss the contributions of the 0-1-jump to the changes in the Oj which we divide in three classes. 
The first class contains the intra-clique contributions which are all contributions coming from sites in the clique where 
the jump takes place (i.e. the central clique in Fig. |3J). For the (p — 1) empty 0-sites in this clique the number of 
occupied neighbouring sites remains unaltered through the jump. A non-zero intra-clique contribution is due to the 
occupied 0-site. After the jump, this site is empty with exactly one occupied neighbouring site which means that the 
number of sites corresponding to must be increased by one. The second intra-clique contribution stems from the 
ending site which has one occupied neighbour before the jump but becomes an occupied site after the jump. This 
means we must decrease the number of vacancies corresponding to by one. 

The second class of contributions will be termed as backward contributions. They stem from the neighbouring 
sites which the particle leaves behind during the jump (labelled with b in Fig. |3J). The backward 0-sites may have 
j = 0, . . . , k occupied neighbours in addition to the jumping particle. The corresponding sites are not shown in Fig. [31 



only the edges are drawn as dashed lines. We will have a positive contribution to from the sites having j occupied 



neighbours in addition to the jumping particle. The same contribution with a negative sign must be attributed to 
o~j+i- In order to quantify the contribution, we need to calculate the expected number of backward sites with j extra 

neighbours. For the calculation we define the partition function si^L which refers to a 0-rooted tree with occupied 
root having exactly Uj neighbouring 0-sites with j extra occupied neighbours. The quantity is given by 



( kll> " ' , ^^e^^ I, H[^(e^e-''r > f)" (^e^e-^E^^ . (27) 



V; 



1=0 ' x /=0 



The average number of backward 0-sites with j additional occupied neighbours is then given by 

fc(p-l) _(0) J + l (0) 

y „ n^ = fc (p_i) mi j+\ . (28) 

Z-~l 1 „(0)' ^ ' x^k 1+1(0) v ' 

u j= ^* 1^1=0 k+1 "l+l 

(0) f 

Here Hi serves as a normalisation as it corresponds to the unconstrained scenario. The r.h.s. is obtained by 

substitution of E^ and S*% according to Eqs. I|20i(l and 1)27(1 and performing the summation. With this result 
the backward contribution of the 0-sites can be easily expressed (see Appendix 0). The calculation of the backward 
contribution of the 1-sites follows very similar lines yielding an average number of backward 1-sites with j additional 
occupied neighbouring sites which is 

i±i (l) 

k l k+li+ \ u . (29) 

X^ k ill J 1 ) V ; 

1^1=0 k+1 °l+l 

A third class contains the forward contributions which come from the sites labelled with f in Fig. The f-sites 
having exactly j — 1 occupied neighbouring sites will yield a positive contribution to and a negative contribution 
comes from the f-sites with exactly j occupied neighbouring sites. The average number of f-sites with exactly j 
occupied neighbouring sites can be calculated with the help of the corresponding partition function 



1=0 



(30) 



One obtains the average number as 



kp k+l-j (0) 



— 1 v fe k+i-i -(o) v ; 

v 3 =0 "0 1^1=0 k+1 u l 
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With similar arguments all the contributions can be derived (for the results see Appendix and one ends up with 
a system of 2(k + 2) differential equations which is of the form 



JO) 



7 (1) 



(32) 



with j = 0, . . . , k + 1 and t = 0, 1. 

For given initial conditions p^(t = 0) = p^ and p^\t = 0) = p^ 1 -* we first calculate the initial values a^ X \t = 0) 
which can be performed following similar lines to the calculation of (1 — pCV 1 ))^ ' (see Sec. IIV Al Eq. (J7J). 
(1 — ( o( / 1 ))'P : [ ^ 1 '' is the probability for an empty site without occupied neighbouring sites which is obviously equivalent 
to ctq ^ 1 "*. A more general calculation yields, for j = 0, . . . , k + 1, 



*r>(-o, = (w-)C;0(>-^^y( 



VP 



(0) 



1 - p(o/i) 



fe+l-j 



(33) 



With these initial values, the <7j-approximation is obtained by integrating Eqs. Ij32|l numerically. We add that the 
conditions l|23|) are fulfilled by the initial values and preserved by Eqs. I|32fl . 

Eqs. H32|) can be evaluated asymptotically as this was done analytically for the p-approximation in Sec. IIV Al The 
Jacobian matrix A of the system is a 2(k + 2) x 2(k + 2)-matrix: 



A : 



^-00 Aqi 

Aw An 



where A 7 - 



\ 



\ 



(m) 
'k + 1 



da. 



(m) 
k + 1 

^ ) 



(34) 



The matrix can be calculated with the help of a computer algebra system as well as its evaluation at a stationary 
points of Eqs. IMl'l) can be done. We find that two eigenvalues always vanish which is due to the conditions l)23|) that 
reduce the effective number of obscrvables by two. If the other eigenvalues {A^} are all negative, the stationary point 
is stable and the inverse relaxation time is given by t _1 = — max{Ai}. 

In the next section we show results for relaxation times and for the time-flow of the observables according to the 
p-approximation and ix, -approximation. 



C. Comparison with MC simulations 

As the exact dynamics is not accessible so far, there is no intrinsic way to check the quality of the results obtained 
according to the PAS. We therefore performed MC simulations implementing the dynamics defined in Sec. Ill Bl oil 
large graphs, details of the generation of the graphs were already given in |l7j . 



Equilibration time 

As a first example, we look at the relaxation time for a system with k = 2 and p = 1, i.e. on an ordinary Bethe 
lattice with coordination number three. On this lattice, due to p = 1, the crystallisation transition is of second order. 
The simulations where performed for an average annihilation rate qo — 0.2 and a mobility of q s — 0.8. The parameter 
qi was used to tune the chemical potential p, according to the condition @ of detailed balance. 

Let us shortly discuss the results of the analytical p-approximation for these parameters. For very diluted systems 
(<7i — > 0, jU — > -co), the relaxation time goes to l/go- For increasing p the relaxation time decreases until it reaches a 
local minimum at p eq ~ —0.445, the point where the slowest process changes from being the decay of global density 
fluctuations to the decay of differences between the sublattice densities. The equilibration time starts to grow until it 
diverges as \p — p cr | _1 if we approach the critical point p cr ~ 1.383. It restarts to decrease and at the end saturates 
at 1/T eq (p ~ * oo) = qo + q s /p which equals one for the above choice of parameters. 

The <7j-approximation leads to qualitatively similar, but quantitatively slightly different results: It starts at the 
same value of the equilibration time for p — > — oo, but reaches the local minimum already at p eq ~ —0.556. The 
behaviour close to the phase transition is extremely similar, with a slightly different prefactor to the dominant critical 
divergence. At large chemical potentials, we have l/r eq (p — > oo) = |go + qs = 0.95, i.e. the equilibration time in the 
vicinity of the closest packing is slightly larger than predicted by the simpler approximation. 
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FIG. 4: Inverse relaxation time for qo = 0.2 and q s = 0.8. The dashed line shows the results of the p-approximation, the full 
line those of the cr., -approximation. MC data is represented by the symbols. The right figure enlarges the critical region to 
show the linear vanishing of 1/r. The slopes on both sides of the transition are related by a factor two. 



For comparison, we have plotted these results together with numerical data in Fig. 01 The numerical simulations 
were performed on a large bipartite lattice of N = 5 ■ 10 6 sites, as an initial condition we have chosen one densest 
packing, the data were averaged over 10 independent runs. For extracting the equilibration time, we have plotted 
the logarithm of the difference p^ (t) — Pg tat from the analytically known equilibrium value, as a function of the MC 
time. The function was linearised in a suitable region. The problem hereby is to identify this region: For short times, 
the system is still far from its asymptotic behaviour, for large time the dynamics is dominated by fluctuations. The 
estimated error of this procedure grows up to 10% for the values close to the transition point. 

Fig. ^demonstrates that both approximations show a very good qualitative agreement with the numerical findings. 
The more detailed Cj-approximation gives, however, the better quantitative values which, within the estimated error, 
almost coincide with the MC results. 



Equilibration of the sublattice densities 

As a further comparison we have investigated the relaxation of the sublattice densities for a system with k = 3 
and p = 2. We have performed the numerical experiments on a graph of N — 4.5 • 10 6 vertices, and have chosen the 
dynamical parameters as q± = 1, qo = e~ 2 and q s = 1 — go. These values correspond to a chemical potential fi = 2 
which is situated in the interval where the liquid as well as the normal crystalline phase are locally stable. 

We have chosen two slightly different initial conditions, in both cases the 0-sublattice was empty. The 1-sublattice 
was first initialised with a partial density p^(t = 0) = 0.145 (experiment A), in the second case with p {1 \t = 0) = 
0.160 (experiment B). After a similar evolution in the initial time interval, the systems approach the two different 
locally stable solutions of the model. Whereas the less dense one converges to the liquid solution, the denser one 
crystallises. The local stability of the liquid phase is demonstrated impressively in case A: After an initial increase, 
the global density falls again toward the liquid value. 

In Fig. the numerical results for the global and both sublattice densities are compared to the p- and the <jj- 
approximation. As in the case of the equilibration times, both approximations reflect correctly the qualitative be- 
haviour. The simpler p-approximation, however, shows a rather pronounced difference in the onset of the final 
relaxation toward the stationary point. This is almost completely cured by the more involved cr^-approximation, 
which shows a quantitatively impressing agreement with the MC data over the whole time interval. 

V. DYNAMICS AT THE GLASS-TRANSITION 

The aim of this section is to detect also the glass-transition in the context of PAS, i.e. through its impact on 
time-dependent, dynamical quantities. Apparently, the observables considered so far (essentially particle densities) 
are not suitable to probe the fundamental difference between the liquid and amorphous state. This is clear as, e.g., 
the difference between the liquid solution densities and the 1RSB solution densities is hardly significant at least for 
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FIG. 5: Global density for two systems (k — 3, p = 2) with initial conditions p ' (t = 0) = and p l '{t = 0) = 0.145 (experiment 
A) resp. p^'(t = 0) = 0.160 (experiment B). The asymptotically reached densities are those of the two metastable states of 
the system, the smaller one corresponding to the liquid, the larger one to the crystalline phase. Numerical data was averaged 
over 100 runs. 




FIG. 6: Global and sublattice densities for the two approximations and the MC simulation. For the corresponding initial 
conditions see caption of Fig. 



an amorphous system in proximity to the glass-transition. It seems more promising to focus on the long lifetime of 
arrested configurations in the glassy state which obviously strongly affects the two-time autocorrelation function 

1 N 

= ~y (ti),m(t 2 ) > *i < *2 j (35) 

i=l 

where ?ij(t) € {0, 1} is the occupation number of site i at time t. In the case of a liquid system we expect C to decay 
quickly with increasing <2 — h to the value Cu — p 2 + (1 — p) 2 that corresponds to uncorrelated configurations. 

For a more symmetric approach we slightly modify the definition of C in Eq. 135JI from a two-time autocorrelation 
to inter-system correlation of two copies of the original model. We introduce 



1 N 

C{t) = ^2 (t).nf (*) ' ( 36 ) 



i=l 
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where n\ (t) and nf (t) are the occupation numbers at time t of two systems indexed I and II that are defined on 
the same graph. Apart from their topological structure the systems can be different, in particular we will consider 
the case where both systems are characterised by different dynamical rates go, qi and q s . We recover Def. l|3*5)) for a 
fixed first system with n\ it) = nf (ti), and by setting t = £2- The correlation function C(t) is the quantity we shall 
approximate in the following. 



A. The correlated-p-approximation 

The use of C(t) as a single observable for PAS seems insufficient for the treatment of correlated systems. We 
therefore transport the simplest description of a liquid system which is given by the particle density p to the level 
of two correlated systems / and // which are defined on isomorphic lattices (i.e. every site has a partner site in the 
other system). This can be done by counting the number of isomorphic pairs of sites that are (a) both empty, (b) 
occupied in system / but empty in system //, (c) empty in system / and occupied in system // and (d) occupied in 
both systems. The corresponding densities w.r.t. to the number of isomorphic pairs (i.e. N) are denoted by p ee , p est , 
p* e and p,,. They constitute the minimal set of observables we shall apply PAS to. The correlation function C{t) is 
then obtained via 

C(t)= Pee (t)+p„(t). (37) 

In a first step we introduce the generalised pseudo-equilibrium by defining three conjugate chemical potentials /x e „, 
/i» e and p^. A chemical potential p ee is not required as p ee = 1 — p est — p* e — p** is not an independent observable. 
Analogously to the previous section, we introduce conditional partition functions S ee , H e „, S„ e and S*, which refer 
to isomorphic rooted trees with an isomorphic pair on the root that is given by the index. The iteration of rooted 
trees can be performed as usual with the difference that we must iterate simultaneously in the two systems now. In 
terms of conditional partition functions this reads 



(H. 



, +pZe*Z p - 1 +pZ*eZ P e - 1 +P^lZ 1 +p(p-l)Z e ^ e ~P- 2 f (38i) 

(5^+pS.eSE- 1 )* (38ii) 
(E^+pE^- 1 )" (38iii) 



-(Hi)' 



(38iv) 



We only explain Eq. I38il i.e. the case where the root of the resulting tree is empty in the two systems. This implies 
the following possibilities for each pair of adjacent cliques (content of large brackets): (1) both cliques are empty; (2) 
clique / is empty and clique II carries a particle; (3) as (2) with I «-> //; both cliques carry a particle which are on 
isomorphic sites (4) or on sites that are not isomorphic (5). The contributions must by weighted according to the 
number of realisations: (1) + p- (2) +p - (3) +p - (4) + p(p — 1) • (5). 
For the calculation we introduce intensive quantities (local fields): 



x = —, y = — , z = — . (39) 

^ee ^ee ^ee 

We can obtain equations for the local fields by taking the corresponding ratios of Eqs. (|38|1 . These read 

(40i) 
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V = — — 7—1 (40ii) 

(40iii) 

They still contain the unknown chemical potentials which we eliminate as described in previous calculations via the 
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introduction of physical partition functions (i.e. the substitution k — > k + 1). This yields the equations 

(41i) 



Pe* 
Pee 

P*e 
Pee 

Pee 
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py) k+1 




(1- 


\-px- 


\-py- 


f pz 


+ p(p- 
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(1- 


Ypx- 


\-py- 
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+ p(p- 


l)xy) k+1 



Substituting the chemical potentials in Eqs. 140J1 according to Eqs. 141|) we find 



(41ii) 
(41iii) 

(42i) 

(42ii) 
(42iii) 

This results in a quadratic equation for x when we substitute y and z in Eq. I|42i|) according to Eqs. <j42iifl and 
(|42iii|) . Only one of the two solutions for x is consistent with the constraints that apply to the observables (positivity 
and restriction of the particle density by ^y). We finally end up with a unique solution for x, y and z which is 



Pe* 






x(l 


+ py) 




Pee 


1 - 


Vpx- 


\-py + 


pz + p(p - 


l)xy 


P*e 








+ px) 




Pee 
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\-py + 
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P** 
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Pee 


1 - 


Vpx - 


\-py + 


pz + p(p — 


l)xy 



y/1 + Apafj +2pa-l 

2p(l+/3-pa) 1 1J 
^l + Apafi + 2p(3-l 

2p(l + a-pP) (43U) 

7(1 + px + py + p{p - l)xy) , (43iii) 



where 



Pe* a P*e j P** r a a\ 

a = , (3 = and 7 = . (44) 

Pee PP** Pee PP** pee PP** 

For the derivation of the differential equations for the observables we only discuss four typical cases. The described 
actions always refer to system /: 

• Insertion on an ee-site: The partition function for an ee-site inside the graph is obtained by substituting 
k — ► k + 1 in Eq. I|38il) . For a successful insertion we must require that all neighbouring sites of the ee-site 
in system / are empty which corresponds to the partition function (S^ e + pS e »S^ c T 1 ) fc+1 . The corresponding 
probability is obtained by dividing by the partition function for an unconstrained ee-site: 

(SL+pS^Sf" 1 )^ 1 / 1+px \ k+1 



(3f e +pS e *Hfe 1 +p3» e 3? e 1 +pE^E p ee 1 + p(p - l)S e *3*e£fe 2 ) fe+1 \1 + px + py + pz + p(p — l)xy 

• Insertion on an e*-site: The partition function for an e*-site follows from Eq. I|38ii|) with k — > k + 1. The 
restriction that all neighbouring sites are empty has the partition function e M " (Sf e ) fc+1 leading to a probability 
for a successful insertion which is 

eM£L) fc+1 ( i \ k+1 



eMsL+ps^- 1 )^ 1 \ 1 + py 



• Jump from a *e-site to an e*-site: The partition function for an unconstrained *e-site follows from Eq. (|38iii(l 
with k — > k + 1. We ask for the probability that the ending site is e* (and not ee). The corresponding partition 
function is (S^Sf" 1 ) 3* e leading to the probability 

\^e*^ee J ^*e _ ^ 



("L + P ^>e*^ee ) ^ P X 
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for an ending site which is e*. For a successful jump the ending site in system / must not have any other 
occupied neighbouring sites. The corresponding partition function is e Mc * (S* e Sf~ 1 )(Sf e ) fc where the ending 
e*-site is regarded as the central site (exponential factor) and the first brackets stand for the clique where the 
jump is performed. The unconstrained case corresponds to e M " (S* e Sf ~ 1 )(Sf e +pS* e Sf~ 1 ) fc which leads to the 
following probability for a successful jump: 



eMH.eHST 1 )^ +pZ*eE P ee 1 ) k U + PV ) 



• Jump from a **-site to an ee-site: In this case the ending site is always ee as the isomorphic partner of the 
starting site is also occupied. The general partition function describing the situation is (S*„H§ ~ x ) E ee ' where 
the ending site chosen as the central site and the first brackets correspond again to the clique where the 
jump is performed. The case where the ending site has no other occupied neighbouring sites corresponds to 
(S^Sf" 1 ) (E p ee + S e *Sf~ 1 ) fe leading to a probability for a successful jump which is 

(5»5g- 1 )(5g e +p5 e .5g- 1 ) fc = / l+px 

e _1 ) Ej \l+px+py + pz+p(p- l)xy 




The corresponding expressions for system II can be obtained by exchanging x and y. 

For the construction of the rate-equations for the observables we must check how the different actions of the 
dynamics contribute to the change in the observables. With the probabilities for the actions which are calculated as 
shown above, we finally obtain 



pee = +P*e <7q ~ Pee q[ 
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(45i) 



(45ii) 



(45iii) 
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-p** qo + p e * q{ { r— - — ) 
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(45iv) 



Here q 1 ^ 11 q\^ 11 and q l J 11 denote the probabilities for particle annihilation, creation and jumps in the respective 
system. 

Applying normalisation p ee + p„ + p* e + p** = 1 and the relations p* e + p»» = p 1 , p e » + p»» = p , where p 1 and 
p u are the particle densities in the respective system, as well as Eq. 13 7|) we find that at any time t: 

Pee{t) = \{C{t) + 1 - /(i) - p"(t)) p„(t) = |(1 - C{t) - /(i) + ^(i)) 

(46) 

p* e (t) = C(t) + p'(t) - p"{t)) p»(t) = i(C(t) - 1 + /(t) + ^(t)) . 

These equations allow for a calculation of the initial values for p ee , p e *, p* e and p** from the given initial conditions 
C(t = 0) = C, p*{t = 0) = p 1 and p = 0) = p 1 . The approximation for the time-flow of the observables is obtained 
from Eqs. I|45() through numerical integration. The correlated- p-approximation for C(t) is then obtained via Eq. (|37|l . 
We add that the approximation was derived following the lines of PAS. 



Stability of the uncorrelated configurations 

After the more general calculation for two correlated systems we return to the study of the autocorrelation function 
of a single system. This can be easily achieved by setting q ! Q = q[ = q% = and qtf — qo, q{ ! — qi and q 1 / = q s in 
Eqs. (|45|l . i.e. we freeze system / in its initial configuration which we assure to be identical to that of system // by 
imposing C(t = 0) = 1. 

We mention that, even though C is a dynamical quantity, the subsequent study does not primarily deal with strictly 
speaking non-equilibrium systems. The interest will be shifted to non-equilibrium situations of the correlation between 
systems that are themselves in equilibrium. Consequently, we assume that system II has a constant particle density, 
i.e. p 1 = p n = p = const. Under these premises Eqs. (I46|) become 

Pee (t) = | (C(t) + 1 - 2p) p,. (t) = \ (C{t) — 1 + 2p) 

(47) 

p„(t) = p. e (t) - i(l-C(t)) . 

Using this result, we obtain a closed differential equation C = C(C) for C{t) when we sum up Eqs. I|45i(l and l|45iv|l . 
One can show that Cu = p 2 + (1 — p) 2 , which corresponds to uncorrelated configurations, is a fixed point of the 
differential equation for C(t). In order to check its stability, we calculate 

dC (1 - (p + 1)P)\(VW - 1)P + 1) (gl(l - - (P + + gap) x (48) 



9C 



C=C„ _ p(l-p)^(pp2 + (l-p)2) 



Here p determines the chemical potential p via the solution for the liquid phase and q = qi/e^ is given through 
detailed balance. 

The fraction on the r.h.s. is always positive, the sign of the expression only depends on the last factor. We conclude 
that the fixed point Cu is stable for 

and unstable for p > p g . Here p g coincides with the spin-glass instability calculated in |17| which is the local instability 
of the liquid solution toward a 1RSB, i.e. glassy, solution. In the case of a continuous transition (p = 1), p g is the 
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FIG. 7: The inverse asymptotic relaxation time of C(t) toward Cij from numerical simulations (upper bounds) compared with 
the relaxation time from the correlated-p-approximation. 



transition density from the liquid to the spin-glass phase. The result appears very reasonable, as we expect the system 
to become, physically speaking, very "viscous" when the density approaches p g coming from the liquid phase. This 
should result in a divergence of the relaxation time of C(t) to Cu which is implied in the stability result for the fixed 
point. We shall compare the results for the relaxation time with the values extracted from MC simulations in Sec. lVBl 
In the case of a discontinuous transition (p > 1) the situation is more complex (see 0,^3)- The relaxational 
dynamics is governed by metastable glassy states that appear for densities below the spin-glass instability at p g . The 
result for the stability of C\\ is unphysical then as the divergence of the relaxation time should occur at lower densities 
(depending on initial conditions etc.). This inconsistency is, however, not surprising as we can show by a 1RSB 
treatment of the ensemble given by with the observables p ee , p e *, p* e and p** that, in the discontinuous case 
(p > 1), the vicinity of C\\ is not described correctly by a replica symmetric approach (which is the approach that was 
used for the correlated- p-approximation) when the density is sufficiently close to p g . We shall therefore only consider 
the continuous glass-transition (p = 1), i.e. the ordinary Bethe lattice, in the following. 

B. Comparison with MC simulations 

We have performed MC simulations for a lattice [k = 2, p = 1) with N = 5 ■ 10 6 sites which was realised as a 
regular random graph. A given dynamics (qo, qi,q s ) was performed for a time T/2 in order to equilibrate the system. 
Subsequently the autocorrelation function C(t) was recorded w.r.t. t\ = T/2 for another time interval of T/2. In 
order to guarantee equilibration in the first half of the simulation time, we have to choose T 3> r eq . To achieve this, 
we have fixed T such that the decay of C{t) to C\\ could be observed. The dynamics was performed without particle 
movement, i.e. qo = min(l,e~''), q± = min(l,e''), q s = 0. 

Asymptotic relaxation time 

From the decay curves of C(t) we can extract asymptotic relaxation times through a linearization of the function 
f(t) = ln(C(t) — Cn) near t ~ T. From the simulations we obtain f(t) as a convex function, which means that the 
extracted relaxation times must be considered as lower bounds for the actual relaxation times. Results for the inverse 
asymptotic relaxation time r _1 are shown in Fig. [7] The numerical data are plotted together with the result of the 
correlated- p-approximation, which follows directly from Eq. I|48l) . 

As mentioned before, the approximation correctly predicts a divergence of the equilibration time at p, g ~ 3.34. 
The quantitative agreement is, however, poor when we approach the phase transition, only fare away numerical and 
approximate analytical data are close to each other. The figure indeed indicates that also the critical exponent of the 
equilibration-time divergence is different for the two cases. This suggests that the correlated p-approximation is still 
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missing the asymptotically dominating slow process, i.e. a more detailed and well chosen set of order parameters is 
needed. 



Decay of the autocorrelation function 

In order to complete the picture, we show the decay curves of C(t) for different /i. The time axis has been shifted 
as to map T/2 — > 0. 

The left of Figs. |H1 shows the curves for small /i, i.e. far away from the spin-glass transition point. The qualitative 
and quantitative agreement between numerical and analytical data is quite satisfactory, even if the differences grow 
with larger values of the chemical potential. Note in particular that the approximation results are very good for 
short times, where by construction the equipartition assumption of PAS is valid, and for large times, where again by 
construction the assumption becomes correct. The largest deviation is found for intermediate times. 
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FIG. 8: Decay curves of the autocorrelation function C(t) for various /x. 



The right of the figures shows analogous curves for positive [i. Since we are coming closer to the transition, the 
deviations show up more drastically, and the quantitative agreement for intermediate times becomes quite poor. Still, 
the initial part and the asymptotic value are reproduced correctly by the correlated-p-approximation. 



VI. CONCLUSION AND OUTLOOK 



In this paper, we have discussed the dynamical behaviour of a hard-sphere lattice-gas model on generalised Bethe 
lattices. The paper can be seen as a natural extension of the results presented in a preceding publication |17| which 
was entirely dedicated to the study of equilibrium properties. 

Since a technique allowing for an exact analysis of the non-equilibrium dynamics of finite-connectivity systems is 
still missing, we have applied a projective approximation scheme. Within this scheme, the dynamics is approximated 
by the time evolution of a finite number of global observables. Since the exact equations for their time evolution do not 
close, we have to use an approximate closure scheme. We have argued that the best and most natural closure scheme 
based only on the considered variables is the following: The true non-equilibrium state of the system is approximated 
by a pseudo-equilibrium distribution which is flat over all microscopic configurations being consistent with the values 
of all considered observables. The quality of this approximation depends crucially on the selection of the observable 
set: The latter has to include for sure a set of observables which are able to describe correctly the equilibrium state 
and to identify different thermodynamic phases. In general, a larger observable set gives also quantitatively better 
results, but the complexity of the underlying calculations increases considerably. 
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As a first application, we have studied the dynamics close to the liquid-crystal transition, both in the case of a second- 
and a first-order transition. Already the minimal description shows a very good qualitative and reasonable quantitative 
agreement with Monte-Carlo simulations, a more sophisticated level of approximation is hardly distinguishable from 
the numerical findings. 

As a second application, we have also studied the dynamics close to the continuous spin-glass transition for the 
model on ordinary Bethe lattices. The chosen minimal approach is able to reproduce correctly the phase transition 
point, where the relaxation time diverges, but the quantitative agreement between numerical and analytical data for 
the dynamical behaviour close to the phase transition is quite poor. This signals that the selected observable set is not 
yet able to include the slowest relaxation process, and richer observable sets have to be used. We have in fact tried to 
do so, but the complexity of the resulting calculations was too large to arrive at a satisfactory result. Concerning the 
approximation scheme close to the spin glass transition, one probably should also go to observable sets which allow 
for the inclusion of explicitly replica symmetry breaking fluctuations. It may well be that the presented approach is 
able to reflect the longitudinal instability of the liquid solution within the replica-symmetric framework, and not the 
important replicon fluctuations breaking the replica symmetry [3| ■ Further research is necessary in this direction. 

Another open question concerns the divergence of the equilibration time close to a discontinuous spin-glass transition 
as occurring on the generalised Bethe lattices. The latter is of fundamental interest in the context of the structural 
glass transition, but cannot be seen using PAS for the simple observable sets considered so far. This would also 
be of importance for the applicability of the method to kinetically constrained models |3Ct l3l| , which have a trivial 
thermodynamics, but show a dynamical arrest due to kinetic constraints. 

Acknowledgement: We are grateful to G. Semerjian and A. Zippelius for many interesting discussions. HHG 
acknowledges also the hospitality of the ISI Foundation in Turin, where parts of this research were worked out. 



APPENDIX A: CONTRIBUTIONS TO THE CHANGES IN a 3 

We introduce the following functions of j: 

P r\j) = (riSlS 717 fOT -? =o <---< fc (ad 

I otherwise 

!J+1 g-C/l) 
te ' = —* (A2) 
otherwise 

Pq°^(j) gives the probability that a 0/1-site, which is contained in an empty clique, has j occupied neighbouring 

sites. P$ 0/1) (j) gives the probability that a 0/1-site, which is empty and which is contained in a clique that carries 
exactly one particle, has j more occupied neighbouring sites (in other adjacent cliques). The values are exact in the 
Scj-ensemble. 

The actions of the dynamics are listed in the following together with their probabilities to occur: 



ID 


Action 


Probability 


KO 
JO 

m 
n 

JOO 
JOl 
JIO 


Removal of a particle from the 0-lattice 
Insertion of a particle to the 0-lattice 
Removal of a particle from the 1-lattice 
Insertion of a particle to the 1-lattice 
Jump of a particle from a 0-site to a 0-site 
Jump of a particle from a 0-site to a 1-site 
Jump of a particle from a 1-site to a 0-site 
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The following table contains the contributions of the different actions of the dynamics to the changes in a- or 
more precisely to the changes in the corresponding numbers of vacancies. The contributions are averages which are 
exact in the Sc-ensemble apart from the exact ^-contributions. 



ID 


Contribution to af^ 


Contribution to cr^ 


Nature of the contribution 


no 


+ (fc + l)(p-l)P 1 (0) (j) 
- (k + lHo- 1) P^(i - 1) 


+ (fc + l)P 1 (1) (. ? ) 

- ffc + I s ) p (1 Vi - 


direct 
neighbours 
uoiphhonrs 


10 


- S 0,j 

+ (fc+l)( 23 -l)p (0) (j-l) 
-(fc+l)(p-l)p (0) (j) 


+ (fc + l)P (1) (j-l) 
-(fc + l)P (1) (j) 


direct 
neighbours 
neighbours 


TZl 


+ (fc+l)pp 1 (0) (j) 
-(fc+ljppwtf-l) 


-1- f)rt A 
~ U U,3 


direct 
neighbours 
neighbours 


11 


+ (fc+l) P F (0) (j-l) 

- (k+ l)oPn (0) f?) 


- S 0.j 


direct 
neighbours 
neighbours 


J00 


+ fc(p-l)pW(j) 
-/e(p-l)p(°)(j-l) 
+ k(p-l)P^(j-l) 
-k(p-l)P^(j) 


+ kP^(j) 

+ fcP (1) (j-l) 
-kP^(j) 


backward 
backward 
forward 
forward 


J 01 


+ fc( P -i)p 1 (0) (j) 

-fc(p-l)p(°)(j-l) 

+ fcpp (0) 0'-i) 

-k P P^\j) 


-Sij 

+ fcP 1 (1) ( J ) 
-fcP«(j-l) 


intra-clinno 
backward 
backward 
forward 
forward 


J10 


- $1,3 

+ kpP[°\j) 

-kpP[°\j - 1) 

+ k(p-l)P^(j-l) 

-fc(p-l)P (0) (7) 


+ fcP (1) (j-l) 
~kP^(j) 


intra-clique 
backward 
backward 
forward 
forward 



The differential equations for the time-evolution of the observables crj ^ are obtained as follows. Let N^ 1 ^ 

be the number of vacancies corresponding to cr^ 1 ^ ■ The average change of N^ 1 ^ during the Monte-Carlo-timc 

At = l/N be A7vj / 1 ' ) which is obtained by choosing the contributions from the respective column in the above table, 
multiplying them with the probabilities for the corresponding actions of the dynamics and summing up the terms. 
The equations for Aivj ^ become differential equations for cr^ ^ in the thermodynamic limit (N — > oo) according 
to 

(nm ArCV 1 ) Mrt 0/1) WW) . 




Here iV^ / 1 ) is the number of sites belonging to the respective sublattice. 
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